############################################################
### Disorganized Political Violence:                     ###
### A Demonstration Case of Temperature and Insurgency   ###
### Replication Code:                                    ###
### Temperature Response Function and Deviation Results: ### 
### Baghdad Nighttime Analyses                           ###
### Andrew Shaver, Alex Bollfrass                        ###
### Final Update: 07/02/2022                             ###
############################################################

# Note: these results were generated using:
# 1) MacBook Pro (15-inch, 2016); Processor: 2.7 GHz Quad-Core Intel Core i7; 
# Memory: 16 GB 2133 MHz LPDDR3
# macOS Catalina
# 2) R 3.6.3 GUI 1.70 El Capitan build (7735)
# 3) RStudio Version 1.3.1093
# (Please see below for specific package versions.)

# Temperature Response and Deviation: Baghdad Nighttime Analyses.

# Load required packages:
library(arm)
library(MASS)
library(lubridate)
library(DataCombine)
library(broom)
library(data.table)
library(radiant.model)
library(lfe)
library(ggplot2)
library(ggpubr)
library(stargazer)
library(locfit)

packageVersion("arm")
# ‘1.11.2’
packageVersion("MASS")
# ‘7.3.51.5’
packageVersion("lubridate")
# ‘1.7.8’
packageVersion("DataCombine")
# ‘0.2.21’
packageVersion("broom")
# ‘0.7.6’
packageVersion("data.table")
# ‘1.12.8’
packageVersion("radiant.model")
# ‘1.3.15’
packageVersion("lfe")
# ‘2.8.5’
packageVersion("ggplot2")
# ‘3.3.5’
packageVersion("ggpubr")
# ‘0.4.0’
packageVersion("stargazer")
# ‘5.2.2’
packageVersion("locfit")
# ‘1.5.9.4’

# Load hour-level datasets:
# Iraq:
iq_sig2_baghdad <- read.csv("~/Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/baghdad_hour_level_data.csv", header = T)
iq_sig2_basrah <- read.csv("~/Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/basrah_hour_level_data.csv", header = T)

# Process data:
# Process dates:
iq_sig2_baghdad$date <- ymd(iq_sig2_baghdad$date)
iq_sig2_basrah$date <- ymd(iq_sig2_basrah$date)

# Create combined Iraq-hour dataset:
# First, create city variable:
iq_sig2_baghdad$city <- "Baghdad"
iq_sig2_basrah$city <- "Basrah"

# For Baghdad time series, data end on 2011-10-31; so remove accordingly:
iq_sig2_baghdad <- iq_sig2_baghdad[iq_sig2_baghdad$date < ymd("2011-10-31"), ]

# Then combine:
iq_hour <- rbind(iq_sig2_baghdad, iq_sig2_basrah)

# Next, aggregate nighttime violence to the day:
iq_hour_nighttime <- iq_hour[iq_hour$times %in% 23 | iq_hour$times %in% 0:6, ] %>%
  group_by(date, city) %>%
  summarise(df = sum(df), idf = sum(idf), ied = sum(ied), 
            max = max(TEMP), temp = mean(TEMP), 
            spd = mean(SPD), vsb = mean(VSB), dewp = mean(DEWP), 
            skc = names(which.max(table(SKC))))

# Create squared and cubic terms:
iq_hour_nighttime$temp2 <- iq_hour_nighttime$temp^2
iq_hour_nighttime$temp3 <- iq_hour_nighttime$temp^3

iq_hour_nighttime$max2 <- iq_hour_nighttime$max^2
iq_hour_nighttime$max3 <- iq_hour_nighttime$max^3

# Create week variable:
iq_hour_nighttime$week <- floor_date(iq_hour_nighttime$date, "week")
iq_hour_nighttime$month <- floor_date(iq_hour_nighttime$date, "month")

iq_hour_nighttime <- as.data.frame(iq_hour_nighttime)

# Next, construct temperature deviations measure:
# Write function:
deviation <- function(x){
  admin_list <- list()
  library(locfit)
  library(data.table)
  for(i in 1:length(unique(panel[, x]))){
    temp <- panel[panel[, x] %in% unique(panel[, x])[i], ]
    temp$nday <- 1:nrow(temp)
    t <- locfit(max ~ lp(nday, nn = .05), data=temp)
    temp$max_hat <- predict(t, newdata = temp$nday)
    temp$max_dev <- temp$max - temp$max_hat
    admin_list[[i]] <- temp
    rm(temp)
  }
  admin_list <- rbindlist(admin_list)
  return(admin_list)
}
panel <- iq_hour_nighttime
iq_hour_nighttime <- deviation("city")
iq_hour_nighttime <- as.data.frame(iq_hour_nighttime)

deviation <- function(x){
  admin_list <- list()
  library(locfit)
  library(data.table)
  for(i in 1:length(unique(panel[, x]))){
    temp <- panel[panel[, x] %in% unique(panel[, x])[i], ]
    temp$nday <- 1:nrow(temp)
    t <- locfit(temp ~ lp(nday, nn = .05), data=temp)
    temp$temp_hat <- predict(t, newdata = temp$nday)
    temp$temp_dev <- temp$temp - temp$temp_hat
    admin_list[[i]] <- temp
    rm(temp)
  }
  admin_list <- rbindlist(admin_list)
  return(admin_list)
}
panel <- iq_hour_nighttime
iq_hour_nighttime <- deviation("city")
iq_hour_nighttime <- as.data.frame(iq_hour_nighttime)

# Create lags:
library(plm)
iq_hour_nighttime <- pdata.frame(iq_hour_nighttime, index = c("city", "date"))
for(i in 1:7){
  iq_hour_nighttime[, ncol(iq_hour_nighttime)+1] <-  plm::lag(iq_hour_nighttime$max, i)
  colnames(iq_hour_nighttime)[ncol(iq_hour_nighttime)] <- paste("max_", i, sep ="")
}
for(i in 1:7){
  iq_hour_nighttime[, ncol(iq_hour_nighttime)+1] <-  plm::lag(iq_hour_nighttime$temp, i)
  colnames(iq_hour_nighttime)[ncol(iq_hour_nighttime)] <- paste("temp_", i, sep ="")
}

for(i in 1:7){
  iq_hour_nighttime[, ncol(iq_hour_nighttime)+1] <-  plm::lag(iq_hour_nighttime$max_dev, i)
  colnames(iq_hour_nighttime)[ncol(iq_hour_nighttime)] <- paste("max_dev_", i, sep ="")
}
for(i in 1:7){
  iq_hour_nighttime[, ncol(iq_hour_nighttime)+1] <-  plm::lag(iq_hour_nighttime$temp_dev, i)
  colnames(iq_hour_nighttime)[ncol(iq_hour_nighttime)] <- paste("temp_dev_", i, sep ="")
}

for(i in 1:7){
  iq_hour_nighttime[, ncol(iq_hour_nighttime)+1] <-  plm::lag(iq_hour_nighttime$df, i)
  colnames(iq_hour_nighttime)[ncol(iq_hour_nighttime)] <- paste("df_", i, sep ="")
}
for(i in 1:7){
  iq_hour_nighttime[, ncol(iq_hour_nighttime)+1] <-  plm::lag(iq_hour_nighttime$idf, i)
  colnames(iq_hour_nighttime)[ncol(iq_hour_nighttime)] <- paste("idf_", i, sep ="")
}
for(i in 1:7){
  iq_hour_nighttime[, ncol(iq_hour_nighttime)+1] <-  plm::lag(iq_hour_nighttime$ied, i)
  colnames(iq_hour_nighttime)[ncol(iq_hour_nighttime)] <- paste("ied_", i, sep ="")
}

for(i in 1:7){
  iq_hour_nighttime[, ncol(iq_hour_nighttime)+1] <-  plm::lag(iq_hour_nighttime$spd, i)
  colnames(iq_hour_nighttime)[ncol(iq_hour_nighttime)] <- paste("spd_", i, sep ="")
}
for(i in 1:7){
  iq_hour_nighttime[, ncol(iq_hour_nighttime)+1] <-  plm::lag(iq_hour_nighttime$vsb, i)
  colnames(iq_hour_nighttime)[ncol(iq_hour_nighttime)] <- paste("vsb_", i, sep ="")
}
for(i in 1:7){
  iq_hour_nighttime[, ncol(iq_hour_nighttime)+1] <-  plm::lag(iq_hour_nighttime$dewp, i)
  colnames(iq_hour_nighttime)[ncol(iq_hour_nighttime)] <- paste("dewp_", i, sep ="")
}

iq_hour_nighttime <- as.data.frame(iq_hour_nighttime)
iq_hour_nighttime$df <- as.numeric(iq_hour_nighttime$df)
iq_hour_nighttime$temp <- as.numeric(iq_hour_nighttime$temp)
iq_hour_nighttime$max <- as.numeric(iq_hour_nighttime$max)

iq_hour_nighttime <- iq_hour_nighttime[iq_hour_nighttime$city %in% "Baghdad", ]

# Direct Fire:
IQ.w.u.l <- felm(log(df+1) ~ max + spd + vsb + dewp + idf + ied +
                   max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                   df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   skc + week  | 0 | 0, data = iq_hour_nighttime)
IQ.w.u.q <- felm(log(df+1) ~ max + max2 + spd + vsb + dewp + idf + ied +
                   max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                   df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   skc + week  | 0 | 0, data = iq_hour_nighttime)
IQ.w.u.c <- felm(log(df+1) ~ max + max2 + max3 + spd + vsb + dewp + idf + ied +
                   max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                   df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   skc + week  | 0 | 0, data = iq_hour_nighttime)

# IED:
IQ.w.i.l <- felm(log(ied+1) ~ max + spd + vsb + dewp + idf + df +
                   max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                   df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   skc + week  | 0 | 0, data = iq_hour_nighttime)
IQ.w.i.q <- felm(log(ied+1) ~ max + max2 + spd + vsb + dewp + idf + df +
                   max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                   df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   skc + week  | 0 | 0, data = iq_hour_nighttime)
IQ.w.i.c <- felm(log(ied+1) ~ max + max2 + max3 + spd + vsb + dewp + idf + df +
                   max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                   df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   skc + week  | 0 | 0, data = iq_hour_nighttime)

# Direct Fire:
IQ.w.u.l.qp <- bayesglm(df ~ max + spd + vsb + dewp + idf + ied +
                          max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                          df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(skc) + as.factor(week), data = iq_hour_nighttime, family = "quasipoisson")
IQ.w.u.q.qp <- bayesglm(df ~ max + max2 + spd + vsb + dewp + idf + ied +
                          max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                          df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(skc) + as.factor(week), data = iq_hour_nighttime, family = "quasipoisson")
IQ.w.u.c.qp <- bayesglm(df ~ max + max2 + max3 + spd + vsb + dewp + idf + ied +
                          max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                          df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(skc) + as.factor(week), data = iq_hour_nighttime, family = "quasipoisson")

# IED:
IQ.w.i.l.qp <- bayesglm(ied ~ temp + spd + vsb + dewp + idf + df +
                          max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                          df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(skc) + as.factor(week), data = iq_hour_nighttime, family = "quasipoisson")
IQ.w.i.q.qp <- bayesglm(ied ~ temp + temp2 + spd + vsb + dewp + idf + df +
                          max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                          df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(skc) + as.factor(week), data = iq_hour_nighttime, family = "quasipoisson")
IQ.w.i.c.qp <- bayesglm(ied ~ temp + temp2 + temp3 + spd + vsb + dewp + idf + df +
                          max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                          df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(skc) + as.factor(week), data = iq_hour_nighttime, family = "quasipoisson")

# First, construct a function for generating counts and creating output that can be plotted.
plot_function <- function(reg_output, X, dataset, type){
  set.seed(1)
  beta.sim <- mvrnorm(100, mu = coef(reg_output), Sigma = vcov(reg_output))
  paste <- model.matrix(reg_output)
  names1 <- colnames(paste)
  beta.sim <- as.data.frame(beta.sim)
  beta.sim <- beta.sim[, c(names1)]
  
  reps.df1 <- matrix(NA, length(unique(round(dataset$temp))), 100)
  for(i in 1:length(unique(round(dataset$temp)))){
    if(X %in% "linear"){
      paste[, 2] <- i + (min(round(dataset$temp))-1)
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta)
      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(X %in% "quadratic"){
      
      paste[, 2] <- i + (min(round(dataset$temp))-1)
      paste[, 3] <- (i + (min(round(dataset$temp))-1))^2
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta)

      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(X %in% "cubic"){
      paste[, 2] <- i + (min(round(dataset$temp))-1)
      paste[, 3] <- (i + (min(round(dataset$temp))-1))^2
      paste[, 4] <- (i + (min(round(dataset$temp))-1))^3
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta)
      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    }
  }
  
  plot_dataset <- as.data.frame(matrix(NA, length(unique(round(dataset$temp))),3))
  colnames(plot_dataset) <- c("means", "cil", "ciu")
  
  plot_dataset$means <- apply(reps.df1, 1 , mean)
  se <- apply(reps.df1, 1 , sd)
  plot_dataset$cil <- plot_dataset$means + qnorm(0.975)*se
  plot_dataset$ciu <- plot_dataset$means - qnorm(0.975)*se
  plot_dataset$temperature <- sort(unique(round(dataset$temp)))
  
  if(X %in% "linear"){
    plot_dataset$model <- "linear"
  } else if(X %in% "quadratic"){
    plot_dataset$model <- "quadratic"
  } else if(X %in% "cubic"){
    plot_dataset$model <- "cubic"
  }
  
  plot_dataset$max <- 0
  plot_dataset[which.max(plot_dataset$means), ]$max <- 1
  plot_dataset$attack_type <- type
  
  return(plot_dataset)
}

# Direct Fire:
p.IQ.w.u.l.qp <- plot_function(IQ.w.u.l.qp, "linear", iq_hour_nighttime, "df")
p.IQ.w.u.q.qp <- plot_function(IQ.w.u.q.qp, "quadratic", iq_hour_nighttime, "df")
p.IQ.w.u.c.qp <- plot_function(IQ.w.u.c.qp, "cubic", iq_hour_nighttime, "df")

# IED:
p.IQ.w.i.l.qp <- plot_function(IQ.w.i.l.qp, "linear", iq_hour_nighttime, "ied")
p.IQ.w.i.q.qp <- plot_function(IQ.w.i.q.qp, "quadratic", iq_hour_nighttime, "ied")
p.IQ.w.i.c.qp <- plot_function(IQ.w.i.c.qp, "cubic", iq_hour_nighttime, "ied")

# Plot:

# Direct Fire:
plot.IQ.m.u.l.qp <- ggplot() + 
  geom_line(data = p.IQ.w.u.l.qp, aes(y = means, x = temperature)) +
  geom_ribbon(data = p.IQ.w.u.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.IQ.w.u.q.qp, aes(y = means, x = temperature), color = "gray4", alpha = 0.5)  +
  geom_ribbon(data = p.IQ.w.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray4") +
  geom_line(data = p.IQ.w.u.c.qp, aes(y = means, x = temperature), color = "gray4", alpha = 0.5)  +
  geom_ribbon(data = p.IQ.w.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray4") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(iq_hour_nighttime$temp)), max(round(iq_hour_nighttime$temp)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(0, 4) +
  ylab("") + 
  xlab("") +
  ggtitle("Direct Fire")

plot.IQ.m.i.l.qp <- ggplot() + 
  geom_line(data = p.IQ.w.i.l.qp, aes(y = means, x = temperature), color = "gray4", alpha = 0.5) +
  geom_ribbon(data = p.IQ.w.i.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.IQ.w.i.q.qp, aes(y = means, x = temperature), color = "gray4", alpha = 0.5) +
  geom_ribbon(data = p.IQ.w.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray4") +
  geom_line(data = p.IQ.w.i.c.qp, aes(y = means, x = temperature), color = "gray4", alpha = 0.5) +
  geom_ribbon(data = p.IQ.w.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray4") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(iq_hour_nighttime$temp)), max(round(iq_hour_nighttime$temp)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(0, 4) +
  ylab("") + 
  xlab("") +
  ggtitle("IED Attacks")

# Temperature deviations:
temperature_range <- 75:105

temperature_deviation_iq_ols <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(df ~ max_dev*threshold + spd + vsb + dewp + idf + ied +
                            max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                            df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            skc + week  | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(df ~ max_dev*threshold + spd + vsb + dewp + idf + ied +
                             max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                             df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             skc + month  | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(df ~ max_dev*threshold + spd + vsb + dewp + idf + ied |
                                 skc + week  | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(df ~ max_dev*threshold + spd + vsb + dewp + idf + ied |
                                  skc + month  | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols_nt <- temperature_deviation_iq_ols(iq_hour_nighttime, 0.95)
iq_90_ols_nt <- temperature_deviation_iq_ols(iq_hour_nighttime, 0.90)

# For plotting Baghdad-Basrah results:
temp_ranges <- unique(iq_95_ols_nt$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges)*2-2), by=2)){
  temp_ranges <- append(temp_ranges, " ", after=i)
}

# Plot:
p_iq_ols <- ggplot() + 
  
  geom_line(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_nt[iq_90_ols_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_nt[iq_90_ols_nt$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_nt[iq_90_ols_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_nt[iq_90_ols_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_nt[iq_90_ols_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_nt[iq_90_ols_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95_ols_nt[iq_95_ols_nt$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90_ols_nt[iq_90_ols_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90_ols_nt[iq_90_ols_nt$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.025, 0.08) + 
  annotate(geom = "text", x = 75:105, y = -0.0225,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# Negative binomial model results:

# Calculate magnitude in terms of percentage deviation from the mean:
new_data_iq <- data.frame(max_dev = c(0, sd(iq_hour_nighttime$max_dev)),
                          threshold = 0,
                          
                          spd = mean(iq_hour_nighttime$spd, na.rm = T), 
                          vsb = mean(iq_hour_nighttime$vsb, na.rm = T), 
                          dewp = mean(iq_hour_nighttime$dewp, na.rm = T), 
                          skc = "CLR",
                          
                          df = mean(iq_hour_nighttime$df),  
                          idf = mean(iq_hour_nighttime$idf), 
                          ied = mean(iq_hour_nighttime$ied),  
                          
                          max_1 = mean(iq_hour_nighttime$max_1, na.rm = T),
                          max_2 = mean(iq_hour_nighttime$max_2, na.rm = T),
                          max_3 = mean(iq_hour_nighttime$max_3, na.rm = T),
                          max_4 = mean(iq_hour_nighttime$max_4, na.rm = T),
                          max_5 = mean(iq_hour_nighttime$max_5, na.rm = T),
                          max_6 = mean(iq_hour_nighttime$max_6, na.rm = T),
                          max_7 = mean(iq_hour_nighttime$max_7, na.rm = T),
                          
                          df_1 = mean(iq_hour_nighttime$df_1, na.rm = T),
                          df_2 = mean(iq_hour_nighttime$df_2, na.rm = T),
                          df_3 = mean(iq_hour_nighttime$df_3, na.rm = T),
                          df_4 = mean(iq_hour_nighttime$df_4, na.rm = T),
                          df_5 = mean(iq_hour_nighttime$df_5, na.rm = T),
                          df_6 = mean(iq_hour_nighttime$df_6, na.rm = T),
                          df_7 = mean(iq_hour_nighttime$df_7, na.rm = T),
                          
                          idf_1 = mean(iq_hour_nighttime$idf_1, na.rm = T),
                          idf_2 = mean(iq_hour_nighttime$idf_2, na.rm = T),
                          idf_3 = mean(iq_hour_nighttime$idf_3, na.rm = T),
                          idf_4 = mean(iq_hour_nighttime$idf_4, na.rm = T),
                          idf_5 = mean(iq_hour_nighttime$idf_5, na.rm = T),
                          idf_6 = mean(iq_hour_nighttime$idf_6, na.rm = T),
                          idf_7 = mean(iq_hour_nighttime$idf_7, na.rm = T),
                          
                          ied_1 = mean(iq_hour_nighttime$ied_1, na.rm = T),
                          ied_2 = mean(iq_hour_nighttime$ied_2, na.rm = T),
                          ied_3 = mean(iq_hour_nighttime$ied_3, na.rm = T),
                          ied_4 = mean(iq_hour_nighttime$ied_4, na.rm = T),
                          ied_5 = mean(iq_hour_nighttime$ied_5, na.rm = T),
                          ied_6 = mean(iq_hour_nighttime$ied_6, na.rm = T),
                          ied_7 = mean(iq_hour_nighttime$ied_7, na.rm = T),
                          
                          week = ymd("2008-05-25"),
                          month = median(iq_hour_nighttime$month),
                          city = "Baghdad")

temperature_deviation_iq_nb <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(df ~ max_dev*threshold + spd + vsb + dewp + idf + ied +
                              max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                              df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(skc) + as.factor(week), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_hour_nighttime$df)
    
    reg_temp_month.nb <- glm(df ~ max_dev*threshold + spd + vsb + dewp + idf + ied +
                               max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                               df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(skc) + as.factor(month) , data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(df ~ max_dev*threshold + spd + vsb + dewp + idf + ied +
                                   as.factor(skc) + as.factor(week), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(df ~ max_dev*threshold + spd + vsb + dewp + idf + ied +
                                    as.factor(skc) + as.factor(month) , data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_nb_nt <- temperature_deviation_iq_nb(iq_hour_nighttime, 0.95)
iq_90_nb_nt <- temperature_deviation_iq_nb(iq_hour_nighttime, 0.90)

# Plot:

p_iq_nb <- ggplot() + 
  
  geom_line(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_nt[iq_90_nb_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_nt[iq_90_nb_nt$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_nt[iq_90_nb_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_nt[iq_90_nb_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_nt[iq_90_nb_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_nt[iq_90_nb_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95_nb_nt[iq_95_nb_nt$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90_nb_nt[iq_90_nb_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90_nb_nt[iq_90_nb_nt$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 75:105, y = -0.0315,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# IED:

temperature_deviation_iq_ied_ols <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(ied ~ max_dev*threshold + spd + vsb + dewp + idf + df +
                            max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                            df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            skc + week  | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(ied ~ max_dev*threshold + spd + vsb + dewp + idf + df +
                             max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                             df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             skc + month  | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(ied ~ max_dev*threshold + spd + vsb + dewp + idf + df |
                                 skc + week  | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(ied ~ max_dev*threshold + spd + vsb + dewp + idf + df |
                                  skc + month  | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols_ied_nt <- temperature_deviation_iq_ied_ols(iq_hour_nighttime, 0.95)
iq_95_ols_ied_nt <- temperature_deviation_iq_ied_ols(iq_hour_nighttime, 0.90)

# For plotting Baghdad-Basrah results:
temp_ranges <- unique(iq_95_ols_ied_nt$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges)*2-2), by=2)){
  temp_ranges <- append(temp_ranges, " ", after=i)
}

# Plot:
p_iq_ied_ols <- ggplot() + 
  
  geom_line(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = iq_95_ols_ied_nt[iq_95_ols_ied_nt$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "firebrick") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.025, 0.08) + 
  annotate(geom = "text", x = 75:105, y = -0.0225,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# Negative binomial model results:

temperature_deviation_iq_ied_nb <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(ied ~ max_dev*threshold + spd + vsb + dewp + idf + df +
                              max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                              df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(skc) + as.factor(week), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_hour_nighttime$df)
    
    reg_temp_month.nb <- glm(ied ~ max_dev*threshold + spd + vsb + dewp + idf + df +
                               max_1 + max_2 + max_3 + max_4 + max_5 + max_6 + max_7 +
                               df_1 + df_2 + df_3 + df_4+ df_5 + df_6 + df_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(skc) + as.factor(month) , data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(ied ~ max_dev*threshold + spd + vsb + dewp + idf + df +
                                   as.factor(skc) + as.factor(week), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(ied ~ max_dev*threshold + spd + vsb + dewp + idf + df +
                                    as.factor(skc) + as.factor(month) , data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "max_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
  
}

iq_95_ied_nb_nt <- temperature_deviation_iq_ied_nb(iq_hour_nighttime, 0.95)
iq_90_ied_nb_nt <- temperature_deviation_iq_ied_nb(iq_hour_nighttime, 0.90)

# Plot:

p_iq_ied_nb <- ggplot() + 
  
  geom_line(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ied_nb_nt[iq_90_ied_nb_nt$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ied_nb_nt[iq_90_ied_nb_nt$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ied_nb_nt[iq_90_ied_nb_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ied_nb_nt[iq_90_ied_nb_nt$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ied_nb_nt[iq_90_ied_nb_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ied_nb_nt[iq_90_ied_nb_nt$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = iq_95_ied_nb_nt[iq_95_ied_nb_nt$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = iq_90_ied_nb_nt[iq_90_ied_nb_nt$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = iq_90_ied_nb_nt[iq_90_ied_nb_nt$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "firebrick") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 75:105, y = -0.0315,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# Finally, combine all results:

p1 <- ggarrange(plot.IQ.m.u.l.qp, plot.IQ.m.i.l.qp)
p1 <- annotate_figure(p1,
                      bottom = text_grob(expression("Temperature (°F)"), 
                                         face = "bold", size = 18),
                      left = text_grob("Expected Count", 
                                       size = 18, rot = 90))
p2 <- ggarrange(p_iq_ols + ylab("OLS") + xlab(""), p_iq_ied_ols + ylab("") + xlab(""), 
                p_iq_nb + ylab("Quasi-Poisson") + xlab(""), p_iq_ied_nb + ylab("") + xlab(""), ncol=2, nrow=2,
                labels = c("", " ", "", ""))
p2 <- annotate_figure(p2,
                      bottom = text_grob(expression(gamma*" (°F)"), 
                                         face = "bold", size = 18),
                      left = text_grob(expression(beta~Delta*"Temperature"), 
                                       size = 18, rot = 90, face = "bold"))

# PRODUCES FIGURE 8 IN PAPER:
ggarrange(p1, p2, ncol = 2)

###
### END.
###
